classdef twoDHist < handle
  % Ad hoc class to encapsulate hist2d and related utilities
  % Derived from handle
  % C-tor sig: Hist = twoDHist(axis, x, y, xbins, ybins, levels),
  %     all args required
  %
  % axis: axis on which to place graph
  % x, y: the data to histogram
  % xbins, ybins: parameters to define the bin boundaries
  % levels: histogram contour-plot levels (or number thereof)
  %
  % Authors: Brandon Sackmann, bsac461@ecy.wa.gov
  %          David Goldsmith dgol461@ecy.wa.gov
  %          Wash. State Dept. of Ecology
  %          et al.?
  % Release date: N/A
  
    properties
        axis; hist; plot;
        x; y; xbins; ybins;
        levels;
    end
    
    methods
        function Hist = twoDHist(axis, x, y, xbins, ybins, levels)
            Hist.axis = axis;
            [~, cols] = size(x);
            if cols ~= 1
                x = x';
            end
            [~, cols] = size(y);
            if cols ~= 1
                y = y';
            end
            n          = length(y);
            mYX        = [y, x];

            vYEdge     = linspace(ybins(1), ybins(2), ybins(3));
            vXEdge     = linspace(xbins(1), xbins(2), xbins(3));
            mHist2d    = Hist.hist2d(mYX, vYEdge, vXEdge);
            maxHist    = max(max(mHist2d));
            contourint = 1:1:maxHist;

            nXBins     = length(vXEdge);
            nYBins     = length(vYEdge);
            vXLabel    = 0.5*(vXEdge(1:(nXBins-1)) + vXEdge(2:nXBins));
            vYLabel    = 0.5*(vYEdge(1:(nYBins-1)) + vYEdge(2:nYBins));

            mHist2d_prob = zeros(size(mHist2d));
            for i = 1:length(contourint)
                tmp = find(mHist2d >= contourint(i));
                mHist2d_prob(tmp)= 1 - (sum(mHist2d(tmp))/n);
            end

            ind_tmp      = find(mHist2d_prob >= 0);
            [X, Y]        = meshgrid(vXLabel, vYLabel);
            [Xsm, Ysm]    = meshgrid(linspace(vXLabel(1), vXLabel(end), 1024), ...
                                     linspace(vYLabel(1), vYLabel(end), 1024));
            mHist2d_prob = griddata(X(ind_tmp), Y(ind_tmp),...
                                    mHist2d_prob(ind_tmp), Xsm, Ysm, 'cubic');
            if ~isempty(find(mHist2d_prob))
                [~, h] = contour(Hist.axis, Xsm, Ysm, mHist2d_prob, 0.05:0.05:0.95, 'fill', 'on');
                CM = get(h,'ContourMatrix');
                CM = CM(:,(CM(1,:) > min(x)) & (CM(1,:) < max(x)));
                xdat = [min(CM(1,:)) max(CM(1,:))];
                ydat = [min(CM(2,:)) max(CM(2,:))];
                set(gca, 'XLim', Hist.fix_lims(get(gca, 'XLim'), xdat));
                set(gca, 'YLim', Hist.fix_lims(get(gca, 'YLim'), ydat));
                caxis([0 1]);
                colormap(jet);
                contour(Hist.axis, Xsm, Ysm, mHist2d_prob, levels, 'color', 'k');
            end
        end
        
    end
    
    methods (Static)
        function mHist2d = hist2d(mYX, vYEdge, vXEdge)
            mHist2d = hist2d(mYX, vYEdge, vXEdge);
        end
        
        function newlims = fix_lims(oldlims, data)
            newlims = oldlims;
            range = oldlims(2) - oldlims(1);
          % fix lower lim, if necessary
            if data(1) - oldlims(1) > 0.1*range
                newlims(1) = data(1) - 0.1*range;
            end
            if oldlims(2) - data(2) > 0.1*range
                newlims(2) = data(2) + 0.1*range;
            end
        end
    end

end

